**Protest Participation Across Europe**

set more off
clear
capture log close
log using euro_prot.log, replace

//Define locals for variables used
local dvs0 = "petition boycott dmnstrt" /*Change below also to allow standalone use of subroutines*/
local dvs = "petit2 boycott2 dmnstrt2"
local iv_c = "pres federal bicam meff newdem" 
local iv_cy = "gini_net rgdpl election enpp un_den"
local iv_i = "age educ income gini_netXincome female married employed union"
local iv_is = regexr("`iv_i'", "gini_netXincome ", "")

use "~/Documents/Projects/Data/ESS/ESS1-4/ess1-4.dta", clear
keep country year s003 s025 dweight `dvs0' `iv_is'
replace income = income-3 		/*Center quintiles*/
replace age=99 if age>99 & !mi(age)
gen age2=age^2
foreach depvar in `dvs0' {
	if regexm("`depvar'","[4]$")==0 {
		rename `depvar' `depvar'2
	}
}
rename petition2 petit2
sort s025

drop if country=="Israel" | country=="Turkey" | country=="Russian Federation" | country=="Ukraine"
replace country="United Kingdom" if country=="Great Britain"
replace country="Slovak Republic" if country=="Slovakia"

egen survey=group(s025)
quietly sum survey
local max=r(max)
label def s025 0 "z"
forvalues i=1/`max' {
	quietly tab year if survey==`i', matrow(y)
	quietly tab s025 if survey==`i', matrow(s)
	quietly levelsof country if survey==`i', local(cc) clean
	local ss = s[1,1] 
	local yy = y[1,1]
	label def s025 `ss' "`cc' `yy'", add
}
label val s025 s025
drop survey

egen c=group(s003)
quietly sum c
local max=r(max)
label def s003 0 "z"
forvalues i=1/`max' {
	quietly tab s003 if c==`i', matrow(s)
	quietly levelsof country if c==`i', local(cc) clean
	local ss = s[1,1] 
	label def s003 `ss' "`cc'", add
}
label val s003 s003
drop c
numlabel s003, add
save ess0.dta, replace

use ess0.dta, clear

*Inequality
sort country year
merge country year using "~/Documents/Projects/Global Inequality/ SWIID v3.1/SWIIDv3_1.dta", ///
	 keep(gini_net) _merge(_m1) nokeep
drop if _m1==1

*GDPpc
sort country year
merge country year using "~/Documents/Projects/Data/GDPpc/pwt70_06032011version/PWT70_cysort.dta", ///
	keep(rgdpl) _merge(_m2) nokeep
replace rgdpl=rgdpl/1000

*Elections and Party Pluralism
sort country year
merge country year using "~/Documents/Projects/Data/Election Data/euro_elections.dta", ///
	keep(election enpp) _merge(_m3) nokeep

*Union Density
sort country year
merge country year using "~/Documents/Projects/Data/Union Density/OECD_Unions.dta",  ///
	_merge(_m4) nokeep

replace year=year-1 if _m4==1
sort country year
merge country year using "~/Documents/Projects/Data/Union Density/OECD_Unions.dta", ///
	_merge(_m4a) nokeep update
replace year=year+1 if _m4==1

sort country year
merge country year using "~/Documents/Projects/Data/Union Density/ICTWSS/ICTWSS_Unions.dta", ///
	_merge(_m4b) nokeep update

replace un_den=21.3 if s025==1002006 /*ICTWSS Database*/
replace un_den=21.3 if s025==1002009
replace un_den=63 if s025==1962006
replace un_den=62.1 if s025==1962008
replace un_den=21.5 if s025==2332004
replace un_den=22 if s025==2332007
replace un_den=22 if s025==2332009
replace un_den=41.6 if s025==7052002
replace un_den=41.3 if s025==7052004
replace un_den=41.3 if s025==7052006
replace un_den=41.3 if s025==7052008

*Institutions
gen pres=(s003==196 | s003==250 |s003==616)

gen federal=(s003==36 | s003==40 | s003==124 | s003==276 | s003==380 | s003==702| s003==756  ///
	| s003==840)

gen bicam=(s003==36 | s003==56 | s003==124 | s003==250 | s003==276 | s003==380 | s003==392 ///
	| s003==528 | s003==616 | s003==702 | s003==724 | s003==756 | s003==840)

*Electoral System
gen meff=.
replace meff=12.5 if s003==40
replace meff=24 if s003==56
replace meff=12.5 if s003==100
replace meff=27.9 if s003==196
replace meff=10 if s003==203
replace meff=25 if s003==208
replace meff=20 if s003==233
replace meff=13.3 if s003==246
replace meff=1 if s003==250
replace meff=10 if s003==276
replace meff=1 if s003==300
replace meff=10 if s003==348
replace meff=4 if s003==372
replace meff=1.3 if s003==380 & year>=2001
replace meff=10 if s003==428
replace meff=15 if s003==442
replace meff=74.6 if s003==528
replace meff=12.5 if s003==578
replace meff=7.5 if s003==616
replace meff=11.3 if s003==620 & year>=1990
replace meff=10 if s003==642
replace meff=10 if s003==703
replace meff=11 if s003==705
replace meff=6.7 if s003==724
replace meff=12.5 if s003==752
replace meff=7.7 if s003==756
replace meff=1 if s003==826

gen lgmeff=log(meff)

*New Democracy
gen newdem=(s003==100 | s003==191 | s003==203 | s003==233 | s003==348 | s003==428 | s003==616  ///
	| s003==642 | s003==703 | s003==705)


sort s025
save ess1.dta, replace

//New hlm commands
local dvs = "petit2 boycott2 dmnstrt2" /*To allow standalone use of subroutine*/
if "`iv_c'"=="" local iv_c = "pres federal bicam meff newdem" 
if "`iv_cy'"=="" local iv_cy = "gini_net rgdpl election enpp un_den"
local iv_i = "age educ income female married employed union"
local iv_is = "age educ income gini_net_X_income female married employed union"

foreach dv in `dvs' {
	use ess1.dta, clear
	capture mkdir `dv'
	quietly cd `dv'
	hlmmkmdm3 using `dv', id2(s025) id3(s003) dv(`dv') l1(`iv_i' dweight) l2(`iv_cy') l3(`iv_c') replace
	hlm hlm3 `dv' int(int(int `iv_c' rand) `iv_cy' rand) age educ  ///
		income(int gini_net rand) female married employed union rand, cmd("`dv'") mdmfile("`dv'.mdm")  ///
		 wgt1(dweight) bernouli store("pa") robust out("`dv'.txt") replace run

	estimates store res_`dv'

	use "`dv'_l1.dta", clear
	sort s003 s025
	merge s003 s025 using "`dv'_l2.dta", nokeep
	drop _merge
	sort s003
	merge s003 using "`dv'_l3.dta", nokeep
	drop _merge
	gen cons=1
	gen gini_net_X_income = gini_net*income
	
	local ivs = "cons `iv_c' `iv_cy' `iv_is'"
	local total_ivs = wordcount("`ivs'")		
	forvalues i = 1/`total_ivs' {
		local t = word("`ivs'", `i')
		local `t'_place = `i'
	}
	
	local stub "b_"
	estsimp2, sims(2500) genname(`stub')
	save "`dv'_1.dta", replace
	
	//Create variables to store median values of all variables	
	use "`dv'_1.dta", clear

	local i=0
	
	*Country vars
	collapse cons `iv_c', by(s003)
	
	foreach a of varlist cons `iv_c' {
		local ++i
		centile `a'
		gen v_`i'=r(c_1)
	}
	
	sort s003
	save "`dv'_3.dta", replace
	
	*Country-year vars
	use "`dv'_1.dta", clear
	
	collapse `iv_cy', by(s025)
	
	foreach a of varlist `iv_cy' {
		local ++i
		centile `a'
		gen v_`i'=r(c_1)
	}
	
	sort s025
	save "`dv'_2.dta", replace
	
	*Individual vars
	use "`dv'_1.dta", clear
	
	sort s003
	merge s003 using "`dv'_3.dta", _merge(_m3)
	sort s025
	merge s025 using "`dv'_2.dta", _merge(_m2)
	drop _m*
	
	foreach a of varlist `iv_is' {
		local ++i
		centile `a'
		gen v_`i'=r(c_1)
	}
	
	save "`dv'_res.dta", replace
	
	//Generate baselines (no gini_net term in base1 and income interaction zeroes out)
	gen base=0
	forvalues i=1/`total_ivs' {
		replace base=base+(v_`i' * `stub'`i')
	}
	
	gen base1=0
	forvalues i=1/`total_ivs' {
		if `i'~=`gini_net_place' {
			replace base1=base1+(v_`i' * `stub'`i')
		}
	}
	
	save "`dv'_res.dta", replace
	erase "`dv'_2.dta" 
	erase "`dv'_3.dta"
	
	
	//table
	capture log close `dv'_results
	log using "`dv'_results.log", replace name("`dv'_results")
	qui sum gini_net
	local min = r(min)
	local max = r(max)
	
	forvalues r = -2/2 {
		local rr = `r'+3
		qui gen est1lo_r`rr' = base1+`min'*`stub'`gini_net_place'+`r'*`stub'`income_place'+`r'*`min'*`stub'`gini_net_X_income_place'
		qui gen est1hi_r`rr' = base1+`max'*`stub'`gini_net_place'+`r'*`stub'`income_place'+`r'*`max'*`stub'`gini_net_X_income_place'
		qui gen pplo_r`rr'=1/(1+exp(-est1lo_r`rr'))
		qui gen pphi_r`rr'=1/(1+exp(-est1hi_r`rr')) 
	
		qui gen diffpp_r`rr'=pphi_r`rr'-pplo_r`rr'
		
		di in white "for inc = " `r'
		sum pplo_r`rr' pphi_r`rr' diffpp_r`rr'
		drop est1lo_r`rr' est1hi_r`rr'
	}
	
	*for comparison
	local j=0
	foreach iv in `ivs' {
		local j=`j'+1
		if `j'~=`gini_net_place' & `j'~=`income_place' & `j'~=`gini_net_X_income_place' {
			qui sum `iv'
			local min = r(min)
			local max = r(max)
			
			qui gen est1lo_`j'=base-v_`j'*`stub'`j'+`min'*`stub'`j'
			qui gen est1hi_`j'=base-v_`j'*`stub'`j'+`max'*`stub'`j'			
			qui sum `iv'
			
			qui gen pplo_`j'=1/(1+exp(-est1lo_`j'))
			qui gen pphi_`j'=1/(1+exp(-est1hi_`j'))
	
			qui gen diffpp_`j'=pphi_`j'-pplo_`j'
			
			di in white "`iv'" ", " `min' " to " `max'
			sum pplo_`j' pphi_`j' diffpp_`j'
			drop est1lo_`j' est1hi_`j'
		}
	}
	log close `dv'_results

	//create data for graph
	*Based on Braumoeller's interaction routine
	qui sum gini_net
	local min=r(min)
	local max=r(max)
	local npts=15
	local inc=(`max'-`min')/(`npts'-1)
	matrix foo2 = 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0  /*One zero for each parameter estimated, plus one*/
	local iter = 0
	while `iter' < =`npts' {
		local giniT = `min'+(`inc'*`iter')
		forvalues i = -2/2 {
			local ii = `i'+3
			qui gen est1`ii'=base1+`giniT'*`stub'`gini_net_place'+`i'*`stub'`income_place'+`giniT'*`i'*`stub'`gini_net_X_income_place'
			qui gen ppest`ii'=1/(1+exp(-est1`ii'))
			qui sum ppest`ii'
			scalar ppe`ii'=r(mean)
			scalar ppl`ii'=ppe`ii'-1.95*r(sd)
			scalar ppu`ii'=ppe`ii'+1.95*r(sd)
		}
		*save values to matrix	
		matrix foo = ppe1, ppl1, ppu1, ppe2, ppl2, ppu2, ppe3, ppl3, ppu3, ppe4, ppl4, ppu4, ppe5, ppl5, ppu5, `giniT'
		matrix foo2 = foo2 \ foo
		drop est* ppest*
		local iter=`iter'+1
	}
	
	*convert matrix to data
	matrix pts=foo2[2..(`npts'+1),1..16]
	svmat pts

	save "`dv'_res.dta", replace
	cd ..
}	

//Table of Coefficients
esttab res_petit2 res_boycott2 res_dmnstrt2 using results_table,  ///
	 not se wide star(* .05) obslast replace

//Predicted Probability Graphs
if "`dvs'"=="" local dvs = "petit2 boycott2 dmnstrt2" /*To allow standalone use of subroutine*/
local lab_petit2=`""Predicted Probability of" "Signing a Petition""'
local lab_boycott2=`""Predicted Probability of" "Boycotting a Product""'
local lab_dmnstrt2=`""Predicted Probability of" "Demonstrating Lawfully""'

//Predicted Probability Graph
foreach depvar in `dvs' {
	cd `depvar'
	use "`depvar'_res.dta", clear
	forvalues i=1(2)5 {
		local m = `i'*3-2
		local l = `i'*3-1
		local u = `i'*3
		twoway connect pts`m' pts`l' pts`u' pts16, lp(solid dash dash) lc(black gs8 gs8) m(i i i) ///
			legend(off) xtitle("") ///
			ytitle("")  ///
			graphregion(lcolor(white) color(white)) nodraw
		graph save ../pp_`depvar'_`i'.gph, replace
	}
	cd ..	
}
graph combine pp_petit2_1.gph pp_boycott2_1.gph pp_dmnstrt2_1.gph, ///
	l1title("Poorest Quintile", size(small) orientation(horizontal) width(17)) ///
	rows(1) ycommon imargin(0 0 0 0) saving(row3.gph, replace) ///
	graphregion(lcolor(white) color(white)) nodraw

graph combine pp_petit2_3.gph pp_boycott2_3.gph pp_dmnstrt2_3.gph, ///
	l1title("Median Quintile", size(small) orientation(horizontal) width(17)) ///
	rows(1) ycommon imargin(0 0 0 0) saving(row2.gph, replace) ///
	graphregion(lcolor(white) color(white)) nodraw

graph combine pp_petit2_5.gph pp_boycott2_5.gph pp_dmnstrt2_5.gph, ///
	l1title("Richest Quintile", size(small) orientation(horizontal) width(17)) ///
	rows(1) ycommon imargin(0 0 0 0) saving(row1.gph, replace)  ///
	graphregion(lcolor(white) color(white)) nodraw
	
graph combine row1.gph row2.gph row3.gph, rows(3) ycommon imargin(0 0 0 0) ///
	graphregion(lcolor(white) color(white) margin(t+15)) plotregion(lcolor(white) color(white)) ///
	xsize(4) b1title("Income Inequality", size(small)) play(protest_titles)  ///
	saving(combo.gph, replace)
	
graph export combo.pdf, replace


//Results Dotplot
use ess1.dta, clear
set more off
local dvs = "petit2 boycott2 dmnstrt2" /*To allow standalone use of subroutine*/
if "`iv_c'"=="" local iv_c = "pres federal bicam meff newdem" 
if "`iv_cy'"=="" local iv_cy = "gini_net rgdpl election enpp un_den"
local iv_i = "age educ income female married employed union"
local iv_is = "age educ income gini_net_X_income female married employed union"

local ivs = "cons `iv_c' `iv_cy' `iv_is'"
local total_ivs = wordcount("`ivs'")		
forvalues i = 1/`total_ivs' {
	local t = word("`ivs'", `i')
	local `t'_place = `i'
}

local lab_ivs = ""
foreach var of varlist `ivs' {
	local lab_`var' : variable label `var'
	local lab_ivs = "`lab_ivs'"+"`lab_`var''"
}
di "`lab_ivs'"

local stub = "b_"

gen ivno = _n in 1/`total_ivs'
keep ivno
drop if ivno==.
save coeffplot_protest.dta, replace

foreach depvar in `dvs' {
	use "`depvar'/`depvar'_res.dta", clear
	quietly gen i_var = ""
	quietly gen ivno = .
	quietly gen b_`depvar' = .
	quietly gen se_`depvar' = .
	centile gini_net
	gen `stub'`income_place'a = `stub'`income_place'+`stub'`gini_net_X_income_place'*r(c_1)
	
	sum gini_net
	replace gini_net_X_income=(gini_net-r(mean))*income
	
	local i = 1
	foreach v in `ivs' {
		local lab_`v' : variable label `v'
		replace i_var = "`lab_`v''" in `i'
		replace ivno = `i'
		sum `v'
		if `r(min)'~=0 & `r(max)'~=1 local sdX2 = 2*`r(sd)'
		else local sdX2 = 1
		if `i'~= `income_place' sum `stub'`i' 
		else sum `stub'`i'a
		
		replace b_`depvar' = `r(mean)'*`sdX2' in `i'
		replace se_`depvar' =`r(sd)'*`sdX2' in `i'	
	
		local i = `i'+1
	}
	
	keep ivno i_var b_`depvar' se_`depvar'
	drop if b_`depvar'==.
	sort ivno
	save "`depvar'/`depvar'_est.dta", replace
	use coeffplot_protest.dta, clear
	sort ivno
	merge ivno using "`depvar'/`depvar'_est.dta", _merge(_`depvar')
	save coeffplot_protest.dta, replace
}

drop if i_var==""
gen order=_n
replace order=-2 if order==`gini_net_place'
replace order=-1 if order==`gini_net_X_income_place'
sort order
drop if i_var=="cons"
keep i_var b* se*
saveold coeffplot_protest.dta, replace

rsource using coeffplot_protest.R, rpath("/usr/bin/r") roptions(`"--vanilla"')

//Messing around with Stata's alternatives for making dotplots
//use coeffplot_protest.dta, clear
//sencode i_var, gen(iv_order)
//reshape long b_ se_, i(i_var) j(dv) string
//egen dv_order = group(dv)
//recode dv_order 1=2 2=3 3=1
//sort dv_order iv_order
//gen ll = b_-1.96*se_
//gen ul = b_+1.96*se_
//
//
//eclplot b_ ll ul iv_order, horizontal ylabel(1(1)18) supby(dv_order, offset(-.25) spaceby(.25))  ///
//	legend(off) rplottype(rspike) estopts(msize(small)) estopts1(mcolor(black))  ///
//	estopts2(mfcolor(gs8)) estopts3(mfcolor(white)) ciopts(lcolor(black) lw(thin)) ///
//	ytitle("") xtitle("") ///
//	graphregion(lcolor(white) color(white)) xline(0, lpattern(dot))
// Things to do: 1. redo workflow to create the data eclplot wants in the first place (long format
// for multiple models, rather than wide), 2. variable labels, 3. vertical zero-marker (an overlaid
// graph), 4. brackets for controls (another overlaid graph?) 


// Correlations
use ess1.dta, clear
collapse petit2-dmnstrt2 gini_net, by(s025)
pwcorr

//Appendix
use ess1.dta, clear
drop if petit2==.
levelsof s003, local(countries)
gen yrs = ""
foreach c of local countries {
	quietly levelsof year if s003==`c', local(years) separate(", ")
	quietly replace yrs = "`years'" if s003==`c'
}
duplicates drop country yrs, force
keep country yrs
sort country
list, noobs clean
save ess_app.dta, replace

set more on
log close
